function [PortWts,allMean,allCovar] = frontier(Universe,Window,Offset,NumPortfolios,ActiveMap,ConSet,NumNonNaN)
%FRONTIER Rolling efficient frontier.
%   [PORTWGTS,ALLMEAN,ALLCOV] = FRONTIER(UNIVERSE,WINDOW,OFFSET,NPORTS,ACTIVEMAP,CONSET,NUMNONNAN)
%   generates a surface of efficient frontiers given data indicating which
%   stocks are active at each date, portfolio constraints and the number of
%   non missing data points needed for each window showing how asset 
%   allocation influences risk and return over time.   
%
%   [PW,AM,AC] = FRONTIER(UNIVERSE,WINDOW,OFFSET,NPORTS)
%   generates a surface of efficient frontiers showing how asset allocation 
%   influences risk and return over time.   
%
%   Inputs:
%
%      UNIVERSE - Portfolio matrix containing the total return data for a
%                 group of securities.  It is an Mx(N+1) matrix where column 1
%                 contains MATLAB date numbers and the remaining columns 
%                 are the total return data for each security.  N is the
%                 number of securities.
%      WINDOW - Number of periods of data to use to calculate each frontier.
%      OFFSET - Increment in number of periods between each frontier.
%      NPORTS - The number of portfolios to calculate on each frontier.
%      
%   Optional inputs:
%
%
%      ACTIVEMAP - An MxN matrix with boolean elements that correspond to
%                  the UNIVERSE where each element indicates if the asset
%                  is part of the UNIVERSE on the corresponding date.  The
%                  default map is an MxN matrix of ones indicating that all
%                  assets are active at all dates.   Note if constraints,
%                  CONSET, other than the default values are used,
%                  the default ACTIVEMAP must be used.
%      CONSET - Portfolio constraints.  Default constraints are generated by 
%               PORTCONS('Default',NumAssets).  This single constraint
%               matrix is applied to each frontier.  Note if an ACTIVEMAP
%               other than the default map, the default constraints must be 
%               used. 
%      NUMNONNAN - Minimum number of non NaN points for each active asset 
%                  in each window of data needed to perform the
%                  optimization.  The default value is WINDOW - N.
%
%   Outputs: 
%  
%     PORTWGTS - A (Number of Curves)x1 cell array where each element is an 
%     NPORTSx(Number of Assets) matrix of weights allocated to each
%     asset.  Number of Assets = length(UNIVERSE). 
%
%     ALLMEAN - A (Number of Curves)x1 cell array where each element is an 
%     1 x (Number of Assets) vector of the expected asset returns used to
%     generate each curve on the surface.
%            
%     ALLCOV - A (Number of Curves)x1 cell array where each element is an 
%     (Number of Assets) x (Number of Assets) vector of the covariance
%     matrix used to generate each curve on the surface.
%
%   See also PORTCONS, PORTOPT.

%       Copyright 1995-2005 The MathWorks, Inc.  

%Number of points on each frontier must be at least 2
if NumPortfolios < 2
  error('finance:frontier:tooFewPortfoliosOnCurve',...
       ['At least two portfolios per frontier must be calculated. \n',...
        'Please specify NPORTS as 2 or greater.'])
end

%Get number of data points between start and end dates.
NumPoints = size(Universe,1);
if NumPoints < Window
  error('finance:frontier:windowGreaterThanNumPoints','Number of observations must be greater than the WINDOW.')
end

%Determine how many curves will be generated given Window and Offset
NumPeriods = ceil((NumPoints - Window + 1)/Offset);

%Get number of assets in portfolio
NumAssets = size(Universe,2)-1;

%Initialize reference security as empty if not given
if nargin < 5 || isempty(ActiveMap)
  ActiveMap = ones(NumPoints,NumAssets);
end

%ActiveMap and Universe must have same dimensions
if ~all(size(ActiveMap) == [NumPoints NumAssets])
  error('finance:frontier:activeMapSizeMismatch',...
        ['ActiveMap must have %d rows and %d columns to match the \n',...
         'number of observations per asset in the specified Universe.'],...
          NumPoints,NumAssets)
end

%End of data window used for first curve
endindex = NumPoints - NumPeriods * Offset + Offset;

%Initialize outputs
X = nan(NumPeriods,NumPortfolios);
Y = nan(NumPeriods,NumPortfolios);
Z = nan(NumPeriods,NumPortfolios);
PortWts = cell(NumPeriods,1);
allMean = cell(NumPeriods,1);
allCovar = cell(NumPeriods,1);

%set up portfolio constraints
if nargin < 6  || isempty(ConSet)
  %No constraints given, use default  
  ConSet = {portcons('Default',NumAssets)};
  ConsInd = ones(1,NumPeriods);
elseif isa(ConSet,'double')
  %Single custom constraint matrix to apply to each curve 
  ConsInd = ones(1,NumPeriods);
  ConSet = {ConSet};
else
  %Each curve uses unique constraint matrix
  ConsInd = ones(1,NumPeriods);
end

%Initialize number of points needed to complete calculation for each window
%if it is not given
if nargin < 7 || isempty(NumNonNaN)
  NumNonNaN = Window;
end

%Need index of start date in first security
sInd = 1;

%Progress bar
h = waitbar(0,'Building efficient frontier...');

try

  %Generate data for rolling frontier
  for period = 1:NumPeriods
    
    % start of Window period
    startindex = endindex - Window + 1;      
    EEnd = Universe(sInd+(startindex-1)+(Window-1),1);
     
    SubMap = ActiveMap(period,:);
    P = find(SubMap);
    SubData = Universe((sInd+(startindex-1):sInd+(startindex-1)+(Window-1)),P+1);
    SubNum = size(SubData,2);
    for i = SubNum:-1:1
      if sum(isnan(SubData(:,i))) > NumNonNaN
        SubMap(i) = 0;
        P(i) = [];
        SubData(:,i) = [];
      end
    end
    SubNum = size(SubData,2);
    SubConSet = ConSet{ConsInd(period)}([1:2,P+2],[P,end]);    
    
    %Estimate mean and covariance of asset returns
    [Mean,Covar] = ecmnmle(SubData);
    Covar = 0.5.*(Covar + Covar');
    
    %calculate portfolios on efficient frontier
    [PortRisk,PortReturn,PW] = portopt(Mean,Covar,NumPortfolios,[],SubConSet);

    %Map weights into matrix based on ActiveMap
    H = zeros(NumPortfolios,NumAssets);
    j = 1:SubNum;
    i = 1:NumPortfolios;
    H(i,P(j)) = PW(i,j);
    
    %Store data to be used for surface visualization
    X(period,i) = year(EEnd) - 1 + (month(EEnd)/12.0);
    Y(period,i) = PortRisk(i);
    Z(period,i) = PortReturn(i);
    
    %Weights, mean and covariance data from points on all curves
    PortWts{period} = H;
    allMean{period} = Mean;
    allCovar{period} = Covar;
    
    %Shift estimation window by Offset
    endindex = endindex + Offset;      
    
    %Update progress display
    waitbar(period/NumPeriods,h,['Building efficient frontier, Period ' num2str(period) ' of ' num2str(NumPeriods)])
  
  end
  
catch exception
    
  close(h)
  error('finance:frontier:unableToDisplayFrontier',exception.message)

end

%Close progress display
close(h)

titlestring = sprintf('Rolling Efficient Frontiers');

%Display surface
surf(X,Y,Z,'FaceColor','interp','EdgeColor','none','FaceLighting','phong');
ylabel('Std. Dev. of Returns');
zlabel('Mean of Returns');
title(titlestring);
camlight right;
